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Previously,  Chang  reported  a  new  high-order  Conservation  Element  Solution  Element 
(CESE)  method  for  solving  nonlinear,  scalar,  hyperbolic  partial  differential  equations  in 
one  dimensional  space.  Bilyeu  et  al.  have  extended  Chang’s  scheme  for  solving  a  one¬ 
dimensional,  coupled  equations  with  an  arbitrary  order  of  accuracy.  In  the  present  paper, 
the  one-dimensional,  high-order  CESE  method  is  extended  for  two-dimensional  unstruc¬ 
tured  meshes.  A  formulation  is  presented  for  solving  the  coupled  equations  with  the  fourth- 
order  in  accuracy.  The  new  high-order  CESE  methods  share  many  favorable  attributes  of 
the  original  second-order  CESE  method,  including:  (i)  the  use  of  compact  mesh  stencil  in¬ 
volving  only  the  immediate  mesh  nodes  surrounding  the  node  where  the  solution  is  sought, 
(ii)  the  CFL  stability  constraint  remains  to  be  <  1,  and  (iii)  superb  shock  capturing  ca¬ 
pability  without  using  an  approximate  Riemann  solver.  To  demonstrate  the  formulation, 
four  test  cases  are  reported,  including  (i)  solution  of  a  two-dimensional,  scalar  convection 
equation,  (ii)  the  solution  of  the  linear  acoustic  equations,  (iii)  the  solution  of  the  Euler 
equations  for  waves  of  small  amplitudes,  and  (iv)  the  solution  of  the  Euler  equations  for 
expanding  shock  waves.  In  all  calculations,  unstructured  triangular  meshes  are  used.  In 
the  first  three  cases,  convergence  tests  show  the  fourth-order  accuracy  of  the  solutions.  In 
the  last  case,  numerical  results  of  the  fourth-order  scheme  are  superior  than  that  obtained 
by  the  second-order  CESE  method. 


I.  Introduction 

In  this  paper,  we  extend  the  one-dimensional,  high-order  CESE  method  for  solving  two-dimensional 
coupled,  nonlinear  hyperbolic  Partial  Differential  Equations  (PDEs)  by  using  unstructured  meshes.  The 
new  formulation  is  currently  fourth-order  in  accuracy,  but  it  can  be  easily  extended  to  a  higher-order  scheme 
by  including  more  terms  in  the  Taylor  series  expansion.  The  tenet  of  the  CESE  method  is  treating  space 
and  time  in  a  unified  manner  in  integrating  the  model  equations  for  flux  conservation.  In  the  CESE  method, 
the  space-time  domain  is  divided  into  non-overlapping  Conservation  Elements  (CEs) ,  over  which  the  space- 
time  integration  is  performed  to  enforce  flux  conservation.  The  integration  is  facilitated  by  using  Solution 
Elements  (SEs),  which  in  general  do  not  coincide  with  the  CEs.  The  unknowns  are  discretized  inside  each  SE 
by  using  a  Taylor  series  expansion.  Aided  by  the  prescribed  discretization  schemes  for  the  unknowns  inside 
each  SE,  space-time  flux  over  each  surface  of  a  CE  can  be  calculated  and  the  overall  flux  conservation  over 
each  CE  can  be  enforced.  The  result  is  an  explicit  formulation  for  updating  the  unknowns  for  time-marching 
calculation.  Special  features  of  the  CESE  method  include:  (i)  The  unknowns  are  discretized  by  using  the 
Taylor  series  expansion  in  both  space  and  time.  The  order  of  the  Taylor  series  is  also  the  order  of  accuracy 
of  the  method,  (ii)  The  method  has  the  most  compact  mesh  stencil,  which  involves  only  the  immediate 
neighboring  mesh  points  of  the  cell  where  the  solution  is  sought,  (iii)  The  method  is  an  explicit  scheme 
in  time-marching  calculation.  The  stability  criterion  is  CFL  <  1.  (iv)  No  approximate  Riemann  solver  is 
used  and  the  scheme  is  simple  and  efficient,  (v)  The  CESE  method  includes  a  suite  of  numerical  algorithms, 
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which  are  built  based  on  the  a  scheme,  a  non-dissipative  scheme.  Extension  of  the  a  scheme  involves  various 
ways  of  adding  artificial  damping  for  maintaining  numerical  stability  and  for  shock  capturing. 

The  mesh  stencil  of  the  new,  fourth-oder  CESE  method  is  identical  to  that  in  the  second-oder  CESE 
method.  For  completeness,  the  treatment  of  unstructured  triangular  meshes  is  repeated  here.  Figure  1 
is  the  projection  of  the  space-time  computational  domain  to  the  two-dimensional  space.  The  solid  lines 
represents  the  triangular  mesh.  Associated  with  each  triangle,  there  are  three  CEs  and  the  dashed  lines 
are  the  boundaries  of  the  CEs.  The  union  of  three  CEs  forms  a  composite  CE,  and  point  3  is  the  centroid 
of  the  CE  associated  with  cell  j.  Symbol  squares  represent  the  vertices  of  the  mesh.  The  open  circles  are 
the  centroids  of  triangles.  The  crosses  are  the  centroids  of  the  composite  CEs.  The  solution  is  stored  at 
the  triangle  centroid.  Figures  2a  and  2b  show  the  CEs  and  SEs  associated  with  the  solution  point  G(j), 
denoted  by  open  a  circle.  These  CEs  are  three-dimensional  space-time  regions,  over  which  flux  conservation 
is  imposed.  In  each  SE,  the  profiles  of  the  unknowns  are  represented  by  a  Taylor  series  expansion.  For 
complete  description  of  the  geometry  of  CE  and  SE,  please  refer  to. 1 

The  remainder  of  the  present  paper  is  organized  as  follows.  Section  II  present  the  numerical  method  for 
the  time  marching  calculation.  This  section  is  broken  up  into  four  sub-sections:  (i)  calculation  of  the  fluxes, 
(ii)  calculation  of  the  temporal  derivatives,  (iii)  calculation  of  the  even-order  derivatives  of  the  unknowns, 
and  (iv)  calculation  of  the  odd-order  derivatives  of  the  unknowns.  In  Section  III,  we  report  numerical  results 
obtained  by  using  the  new  scheme.  Finally,  we  present  our  concluding  remarks  and  list  the  cited  references. 


II.  Numerical  Method 


The  two-dimensional  model  equations  to  be  solved  are  presented  in  the  following  vector-matrix  form: 


<9U  dFx  dF  y 
dt  dx  dy 


(1) 


where  U  =  («i,  n2,  •  ■  • ,  «m)‘,  and  Fx’y  =  (fx ’v ,  ’v, . . . ,  and  there  are  m  coupled  equations  in  Eq.  (1). 

Each  equation  in  Eq.  (1)  can  be  recast  into  a  divergence  free  condition:  V  ■  h,  =  0,  where  i  =  1, . . . ,  m  and 
h,  =  (/f ,  /f  ,UiY  is  the  space-time  flux  function  of  the  ith  equation.  The  divergence  here  operates  in 
three-dimensional  Euclidean  space.  Aided  by  the  Gauss  theorem,  we  have 


hi(x,y,t)  ■  ds  =  0.  (2) 

The  CESE  method  is  a  numerical  method  to  integrate  Eq.  (2)  while  treating  space  and  time  in  a  unified 
manner.  The  integration  is  facilitated  by  the  discretized  conserved  variables  iq  and  the  fluxes  ff’v  with 
i  =  1, . . . ,  to.  First,  the  space-time  domain  is  divided  into  non-overlapping  Conservation  Element  (CE).  Over 
each  CE,  Eq.  (2)  is  imposed.  The  actual  integration  is  facilitated  by  the  definition  of  Solution  Elements 
(SE),  in  which  m  and  f^’v  are  represented  by  Taylor  series  expansion.  In  deriving  the  formulation  for 
integrating  Eq.  (2),  we  make  the  following  assumptions:  (i)  Inside  an  SE,  the  unknowns  and  flux  functions 
are  discretized  and  they  are  represented  by  using  a  Taylor  series,  (ii)  The  desired  order  of  convergences 
is  even.  In  this  paper,  the  fourth-order  method  will  be  developed,  (iii)  Each  even-order  derivative  has  an 
odd-order  derivative  with  respect  to  each  spatial  direction,  (iv)  The  discretized  variables  obey  the  alternative 
rule  of  differentiation,  e.g.,  {d2Ui/ dxdy)  =  {d2Ui/ dydx)  ,  where  the  superscript  *  denotes  the  discretized 
variables. 

II. A.  The  Taylor  Series  Expansion 

To  proceed,  the  conserved  variables,  u\,  m2,  •  •  • ,  um,  and  the  fluxes,  fi’y,  /2  ’v ,  .  ■  . ,  f%y  are  discretized  by 
using  the  Taylor  series  expansion  in  space  and  time.  Shown  below  is  the  third-order  Taylor  expansion  in 
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(3) 


time  and  a  two-dimensional  space  of  a  conserved  variable  m,  where  1  <  i  <  m, 

Ui(x,  y,  t)  =Ui(xj,yj,  tn)  +  (m)x Ax  +  ( m)vAy  +  (ui)tAt+ 

-j  [(ui)xx Ax2  +  2(ui)xyAxAy  +  (■ Ui)yyAy 2  +  2(ui)xtAxAt  +  2{ui)ytAyAt  +  (m)ttAt2]  + 

^  { {^i)xxxAx  -f*  {ui'jyyyAy  T  \(up  xxy  T  2 (ui')xyx\  Ax  Ay  T  [2 (yj>i)xyy  [vji)yyx\  Ay  Ax-\- 
3(ui)xxtAx2  At  +  3(ui)xttAxAt 2  +  3{ui)yytAy2  At  +  3(ui)yttAyAt2+ 

6(ui)xytAxAyAt  +  ( Ui)tttAt 3}  . 

where  the  anchor  point  of  the  Taylor  expansion  is  ( Xj ,  yj,  tn),  and  the  distance  is  denoted  by  Aa:  =  x  —  Xj, 
Ay  =  y  —  yj,  and  At  =  t  —  tn.  Thus  all  derivatives  are  evaluated  at  ( Xj,yj,tn ).  The  above  Taylor  series  can 
be  written  in  a  more  concise  form: 


N  N-aN-a-b 


Ui(x,y,t)  =  ^2^2  E 


Qa+b+c 


1 


-AxaAybAtc 


(4) 


n  .  dxadybdtc  a\b\c\ 

a— 0  6=0  c=0  y 

where  N  is  equal  to  the  desired  order  of  accuracy  of  the  method  minus  1.  For  the  fourth-order  CESE 
method,  N  =  3.  Again,  all  derivatives  on  the  right  hand  side  of  Eq.  (4)  are  evaluated  at  ( Xj,yj,tn ). 
Moreover,  in  Eq.  (4),  we  invoked  assumptions  (iii)  and  (iv),  which  lead  to  d2u/dxdy  =  d2u/dydx  and 
d3u/dx2dy  =  ( uxxy  +  2uxyx)/3 ,  d3u/dxdy2  =  (2 uxyy  +  uyyx)/3  and  so  on.  Similarly,  the  Taylor  expansion 
of  a  derivative  of  the  unknowns  can  be  written  as 

A  A— a A—a—b 


dcm 


dx1 dyJdtK 


{x,y,t)  =  22J2  E 


dBi 


1 


dxI+adyJ+bdtK+c  a\b\c\ 


AxaAybAtc, 


(5) 


a— 0  6=0  c— 0 

where  C  =  I  +  J  +  K,A  =  N  —  C,  and  B  =  C  +  a  +  b  +  c.  We  remark  that  Eq.  (4)  is  a  special  case  of  Eq.  (5) 
with  A  =  N.  Similarly,  the  Taylor  expansion  of  fluxes  can  be  expressed  as 


dcf \ 


x,y 


A  A— a A—a—b 


dx1 dyJ dtK 


(x,y,t)  =  22Y,  E 


dBf2y 


a=0  6=0  c=0 


dxI+adyJ+bdtK+c  a\b\c\ 


AxaAybAtc, 


(6) 


II. B.  The  Chain  Rule 

To  proceed,  we  show  the  relationship  between  the  derivatives  of  the  fluxes  f*'v  and  the  derivatives  of  the 
conserved  variables  Ui,  where  i  =  1, . . .  ,m.  Since  f['v  are  known  functions  of  m,  aided  by  the  chain  rule, 
the  derivatives  of  f-',v  can  be  related  to  the  derivatives  of  m.  For  the  first-order  derivatives,  we  have 


df-’v  3/r v  dm 

94-1  A  dui  dVi 


z=i 


where  4/i  =  x,y,t.  For  the  second-order  derivatives,  we  have 


d2f , 


x,y 


d^id^/2 


=  E 


dff'v  d 2 


m 


dm  94> \d'$>2 


EE 


d2f*’v  dm  dun 
dmdun  chFi  <94/2 


(7) 


(8) 


Z=1  ”  ““  l—l  n—  1 

where  (4/i,  4^2)  =  (x,  x),  (y,  y),  ( t ,  t),  (x,  y),  (x,  t ),  and  (y,  t).  For  the  third-order  derivatives,  we  have 


d3f ■ 


df? 


<9  4/ 1  <9  4/2  54/ 3 


dm  94,i94'294>3 


=E 

z=i 

m  m 

EE 

1  =  1  72=1 

772  772  772 

EEE 

Z=1  72=1  p=  1 


d3u. 


d2m  dun  +  d2m  dun 


22 112 

dmdun  V94'i94,2  34,3  '  34/i94/3  94/2 

d3f2v  dm  dun  dup 
dmdundup  (94/ 1  <94>2  94/ 3 


d2m  dun 
(94,294,3  94>] 


+ 


(9) 


where  (4-i,  4-2,  4'3)  =  (x,  x,  x),  {y,  y,  y),  (t,  t,  t),  (x,  y,  t),  (x,  x,  y),  (y,  y,  t),  (x,  x,  t),  (2,  y,  y),  (y,  t,  t),  and  (2,  t,  t). 
These  equations  are  recursive  and  can  be  further  extended  for  the  high-order  derivatives.  Aided  by  Eqs.  (7-9), 
the  derivatives  of  f*’v  are  related  to  the  derivatives  of  m. 
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II. C.  Additional  Equations 


We  proceed  to  obtain  additional  model  equations  by  applying  spatial  and  temporal  differentiation  to  the 
original  first-order  governing  equations  Eq.  (1): 


d_ 

dx 

d_ 

dy 

d_ 

at 


dft 

df? 

dx 

dy  ’ 

d2f f 

d2/y 

dxdx 

dydx 

d2ff 

d2n 

dxdy 

dydy 

d2f* 

d2f\ 

dxdt 

dydt ' 

(10) 


Equation  (10)  can  be  extended  to  higher-order  derivatives.  In  general,  a  derivative  of  iq  involving  temporal 
differentiation  can  always  be  found  by  using  the  chain  rule,  i.e.,  Eqs.  (7-9),  in  conjunction  with  the  following 
generic  formulation  extended  from  Eq.  (10): 


dcUi  =  dcf f _ dcf! 

dxIdyJdtK  dxI+1dyJdtK~ 1  dxIdyJ+ldtK~1' 

where  C  =  I  +  J  +  K ,  K  =  1, 2, . . . ,  N,  J  =  0, 1, . . . ,  N  —  K,  and  I  =  0,1 , . . . ,  N  —  K  —  J.  Equation  (11)  is  a 
divergence- free  condition  the  derivatives  of  space-time  flux  vector  h, .  Aided  by  the  Gauss  theorem,  we  have 


/d 3pjfer(*.».<)-*  =  0 


(12) 


Equation  (2)  is  a  special  case  of  Eq.  (12).  To  recapitulate,  all  derivatives  of  f-"v  can  be  expressed  by 
derivatives  of  iq,  and  all  temporal  derivatives  of  m  can  be  expressed  in  terms  of  the  spatial  derivatives  of  u-t . 
Therefore,  as  that  in  the  original  second-order  CESE  method,  the  independent  variables  in  the  high-order 
CESE  method  include  only  the  conserved  variables  it,;  and  their  spatial  derivatives.  Table  1  lists  the  primary 
unknowns  of  the  fourth-order  CESE  method. 


Table  1:  Primary  unknowns  of  a  two-dimensional,  fourth-order  CESE  method. 
Even-order  variables  Odd-order  variables 


Ui 

i)x 

(' ui)v 

{^i)xx 

i)xxx 

i)xxy 

i)xy 

(Vji)xyx 

(' ui)Xyy 

{' ui)yy 

( ui)yyx 

(ui)yyy 

II. D.  Space-Time  Flux  Conservation 

In  this  section,  we  integrate  the  original  model  equations,  i.e.,  Eq.  (2)  as  well  as  the  additional  equations 
Eq.  (12)  to  find  all  even-order  derivatives  of  the  conserved  variables.  For  the  solution  at  the  new  time  step 
by  using  the  fourth-order  CESE  method,  we  integrate  Eq.  (2),  i.e.,  the  original  equation,  for  iq,  and  we 
integrate  the  following  special  cases  of  Eq.  (12): 


d2ht 

dxdx 


■  ds  =  0, 


•  ds  =  0, 


ds  =  0 


(13) 


for  the  solutions  of  ( Ui)xx ,  ( Ui)xy ,  and  ( Ui)yy,  respectively.  As  will  be  shown  in  the  following,  the  integration 
results  in  explicit  formulation  for  the  above  unknowns.  Once  all  the  even-order  unknowns  of  a  given  order 
are  calculated,  all  the  odd-order  derivatives  of  the  unknowns  can  be  found  by  using  a  central-difference  like 
procedure.  Central  differencing  results  in  the  solution  of  ( ut)x  and  (ui)y,  and  central  differencing  (ui)xx 
results  in  numerical  solutions  of  ( Ui)xxx  and  ( Ui)xxy ,  et  cetera.  This  procedure  is  consistent  with  that  in  the 
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original  second-order  CESE  method.  The  remainder  of  this  subsection  will  focus  on  the  calculation  of  flux 
conservation. 

The  even  derivatives  are  found  by  conserving  the  flux  through  each  CE.  The  calculation  includes  three 
steps:  (i)  the  flux  through  the  six  side  faces,  (ii)  the  flux  through  the  bottom  face,  and  (iii)  the  flux  through 
the  top  face.  The  flux  through  the  side  and  bottom  faces  are  calculated  by  using  the  known  solutions  at 
the  previous  time  step  at  neighboring  points  of  the  current  solution  point.  Therefore,  the  calculation  for 
fluxes  through  side  and  bottom  surfaces  are  explicit.  The  flux  through  the  top  surface  is  a  function  of  the 
solution  itself.  In  the  following,  the  derivation  is  presented  in  two  subsections  for  flux  through  side  surface 
and  through  bottom/top  surfaces. 


II.D.l.  Flux  Through  Side  Faces 

Figure  3  shows  the  faces  associated  with  the  solution  point  ( x3x,y3x ).  For  a  two-dimensional  triangular 
mesh,  there  are  three  neighboring  triangles  surrounding  the  solution  point  (x3x,y3x).  There  are  six  side 
surfaces  associated  with  the  CE  centered  at  (x3x:y3x).  The  centroids  of  three  neighboring  triangles  are 
denoted  by  (xv,yv),  (xw,y*>'),  and  (xv,  yv).  We  note  that  points  1',  5',  and  7'  are  not  the  solution  points 
of  the  neighboring  CEs.  They  are  the  centroids  of  the  neighboring  triangles.  In  the  following,  we  denote  the 
solution  points  associated  with  points  1',  5',  and  T  by  point  l'x,  5,x,  and  7,x.  Extended  from  each  of  these 
centroids,  a  pair  of  side  surfaces  exist  as  a  part  of  the  surface  of  the  CE  centered  at  (x3x  ,y3x).  On  a  side 
surface,  the  values  of  the  conserved  variables  and  fluxes  can  be  calculated  by  the  Taylor  series  expansion 
from  the  corresponding  solution  points,  i.e.,  (x^x  ,y ^x),  (x3ix  ,y3>x),  or  (xj/x  ,yvx). 

To  proceed,  we  consider  point  (xv,yv)  and  the  two  associated  side  surfaces  spanned  between  (x2',y2') 
and  (£4', ?/4').  We  will  illustrate  the  calculation  of  one  side  surface  spanned  between  (xv,yv)  and  ( X2',y2 ') 
by  using  the  Taylor  series  expansion  from  point  l,x,  which  in  general  does  not  coincide  with  (xv,yv)  unless 
all  triangles  all  equilateral.  We  first  calculate  the  coefficients  of  the  Taylor  series  representing  the  fluxes. 
Then  we  integrate  the  Taylor  series  over  the  side  surface  to  find  the  flux  through  the  surface.  At  a  location 
(x,y)  on  the  side  surface  spanned  by  (xv,yv)  and  (*2', 2/2')  and  within  tn  —  1/2  <  t  <  tn,  the  flux  f?’v  can 
be  expressed  as 

_  V  y  6  da+b+cfty  AxaAybAtc 
1  2-^i  2-^i  dxadybdtc  a\b\c\  ’ 

o=0  b—0  c= 0  y 


where  Ax  =  x  —  Xyx  ,  Ay  =  y  —  yVx  ,  and  At  =  t  — tra_1/2.  The  superscript  *  denotes  the  discretized  variables. 
It  understood  that  the  derivatives  of  the  flux  function  in  the  integrand  are  calculated  at  (xyx  ,2/1'x).  Next 
we  adopt  the  following  notation  to  facilitate  the  integration: 

Ax  =  a(x2'  -  xv )  +  xv  -  xi'x ,  Ay  =  a/2/2'  -  yv)  +  yv  -  2/i'x ,  A t  =  /3^-, 

where  0  <  a  <  1,  and  0  <  /?  <  1.  A  single  spatial  transformation  variable  is  possible  because  the  projection 
of  the  side  face  to  the  spatial  domain  is  a  straight  line  segment. 

The  unit  normal  vector  on  the  side  surface  is  n  =  [(2/2'  —yv),  —  ( X2 >  —  Xv),Q\/ \J (X2'  —  X\> )2  +  (2/2'  —  yv)2 ■ 
The  differential  area  for  the  integration  dR  =  dsdt ,  where  ds  =  da^J (x2>  —  xv)2  +  (2/2'  —  yv)2  and  dt  = 
At/2d(3  Aided  by  these  definitions,  the  total  flux  through  the  side  face  can  be  expressed  as 


(Fi)  1 


•  n  dR 


1'2'21 
1 


/  V  At  .  .At 

\V2'  -y v)~2>-{X2'  -xv)  —  ,0 


dad/3 


=  [[  \ifi)*,(fi r,M*} 

J  J a,/3—0 

=  (2/2'  -  2/1')^  /  /  {fi)*dadP  -  {x2'  -  Xv)^-  [  [  (fV)*dadfi. 

z  0  z  J  J ot,B=0 


(15) 


The  calculation  in  Eq.  (15)  depends  on  available  solutions  at  point  1 x  at  a  previous  time  step.  The  calculation 
of  fluxes  through  the  other  five  side  surface  is  similar  to  Eq.  (15).  As  such,  the  flux  through  all  side  surfaces 
of  the  CE  associated  with  the  solution  point  ( x3x,y3x )  can  be  readily  calculated: 


(Fi)  side  =  (Fi)  l'2'21  +  (Fi)  l'4'41  +  (^i)5'2'25  +  (^i)5'6'65  +  (^i)7'4'47  +  (Fi)re'67- 


(16) 
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II.  D.  2.  Flux  Through  the  Bottom  Surface 

To  proceed,  we  present  the  calculation  of  the  flux  through  the  bottom  surface  of  the  CE  centered  at  point  3X . 
The  bottom  surface  is  at  the  old  time  step  n  —  1/2  and  thus  all  points  are  denoted  by  the  superscript  '.  The 
calculation  of  the  flux  through  the  bottom  surface  of  the  CE  requires  integration  over  three  quadrilaterals. 
Again,  we  will  only  derive  the  calculation  of  the  quadrilateral  associated  with  points  1'  and  l,x  as  shown 
in  Fig.  4.  In  this  figure,  point  1'  is  the  centroid  of  the  neighboring  triangle  and  point  3'  is  the  centroid  of 
central  triangle.  Points  2'  and  4'  are  the  shared  vertices  of  triangles.  To  proceed,  we  split  the  quadrilateral 
into  two  triangles,  i.e. ,  Al'3'4'  and  Al'2'3'. 

To  facilitate  the  integration,  we  first  perform  coordinate  transformation  that  converts  an  arbitrary  triangle 
into  a  right  triangle.  In  the  following  discussions,  we  will  only  focus  on  triangle  Al'3'4'  as  shown  in  Fig.  4.  The 
(x,y)  coordinates  are  transformed  to  (£,?y)  with  (£y,r}y)  =  (0,0),  (£3', 773')  =  (1,0),  and  (fy,r]y)  =  (0,1). 
The  transformation  equation  is 


dx 

dy 


XH  xv  )  (  dt 

Vi  Vn  J  V  drl 


(17) 


Assume  all  entries  of  the  Jacobian  matrix  are  constant  and  we  have  x j  =  (2:3/  —  xy)/(l  —  0)  =  xy  —  xy, 
Xv  =  (xy  -  Xy),  =  (2/3/  -  yy),  and  y n  =  (yy  ~  yy)-  By  integrating  Eq.  (17)  from  point  l'x,  i.e.,  the 
solution  point  associated  with  the  triangle  centered  at  point  1  at  n  —  1/2  time  step,  to  any  other  location 
(x,y)  inside  Al'3'4',  we  have 


(x  —  Xy 

y-w 


(  X3>  -  Xy 

\  2/3'  -  yy 


where  0  <  rj  <  1  and  0  <  £  <  1  —  77,  and  the  Jacobian  determinant 


of  the  coordinate  transformation  is 


d(x,y) 

d&v) 


|(2;3'  -  Xy)(yy  -  yy )  -  {xy  -  Xy)(y3'  ~  yy)\  . 


Aided  by  the  coordinate  transformation,  the  integration  of  the  space-time  flux  h,  over  triangle  <51'3'4'  can 
be  expressed  as 


(Fi) 


1'3'4'  — 


h  ,  •  n  dR 


1'3'4' 


[I  [(/f)*,  (/?)*,  («i)*]-(0,0,-l)Jded»7 
J 

~JJ  (' Ui)*Jd£v 


(18) 


nl-77  A  A-a 

EE 

a=0  6=0 


v  =  /QQ+b 

EEii  Q  a  Q™  )  +  dxu  +  A;ci)a  +  rjyv  +  Ayx)b  Jdfdy. 

r, — n  h — n  ’  ‘  \  y 


yy. 


where  Ax\  =  xy  —  aqx  and  Ayi  =  yy  —  y±x  ■  Equation  (18)  can  be  readily  calculated  because  all  derivatives 
of  Ui  at  point  l'x,  i.e.,  at  the  previous  time  step  tn_1/2,  are  known.  The  fluxes  through  triangle  II  as  well 
as  through  other  parallelepipeds  associated  with  points  5'  and  7'  can  be  calculated  in  a  similar  manner. 
Thus  the  space-time  flux  through  the  whole  bottom  surface  of  the  CE  associated  with  point  3  can  be  readily 
calculated: 


(Fi)  bot  —  (Fi)yyy  +  (F^yw  +  (F))5'6'3'  +  (Fi)  5'3'2'  +  (Fj)  7'3'6'  +  (Fi)ryy.  (19) 

II.  D.  3.  Flux  Through  the  Top  Surface 

The  flux  through  the  top  surface  of  the  CE  associated  with  point  3  can  be  calculated  in  a  like  manner  as 
that  through  the  bottom  surface  of  the  CE.  The  difference  is  that  the  solution  point,  from  which  the  Taylor 
series  expansion  is  expanded,  is  point  3X.  Instead  of  using  three  solution  points  1',  5',  and  7'  as  that  for  the 
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flux  through  the  bottom  surface  of  the  CE,  only  one  solution  point,  i.e.,  point  3X,  is  used  for  flux  through 
the  top  surface: 


(A)  134  = 


hi  •  n  dR 


134 


II  [(AT,  0,1) 

jj  (ui)*Jd£ri 


(20) 


1  /» 1 — Ti  A  A-a 


— ‘  1  fda+bu. 
< 

a— 0  b— 0 


/0  Jo  dxag™b  )  {&i  +  T1xr]  +  Ax3)a(Zyi  +  i1yri  +  Ay3)b  Jd^drj. 


where  AX3  =  x3  —  x3x  and  Ay3  =  y3  —  y3x.  Similar  to  that  for  the  bottom  surface,  the  total  flux  through 
the  top  surface  is 


(-Fi)top  —  {Fi)l34  +  (Fi)  123  +  (F'i)  563  +  {Fi)  532  +  (Fi)  736  +  (-Pi)  743-  (21) 

Equation  (21)  is  an  implicit  function  of  the  conserved  variable  Ui  and  its  spatial  derivatives  at  point  3X.  For 
the  fourth-order  CESE  method,  the  spatial  derivatives  of  concerned  are  up  to  the  third  order.  According 
to  Table  1,  there  are  total  12  even-  and  odd-order  derivatives  of  tq  involved  in  Eqs.  (20)  and  (21).  In  the 
following,  we  illustrate  the  procedure  of  integrating  Eqs.  (20)  and  (21)  in  several  steps.  In  particular,  we  will 
show  that  all  calculation  steps  involved  are  explicit;  there  is  no  need  to  solve  the  function  implicitly,  e.g., 
inverting  a  matrix. 

First,  we  remark  that  point  3X  is  the  centroid  of  the  hexagonal  area  defined  by  points  125674,  which  is  the 
top  surface  of  the  CE  of  concern.  As  such,  the  integration  of  the  conserved  variable  (m) 3x  ,  i.e.,  the  first  term 
in  the  Taylor  series  expansion,  over  the  entire  top  surface  is  simply  (ui)3x  multiplied  by  the  hexagonal  area, 
i.e.,  {ui)3x  Ahex,  where  Ahex  denotes  the  area.  Moreover,  the  integration  of  the  first-order  derivative  terms  in 
the  Taylor  series  expansion,  i.e.,  terms  involving  (ui)x,  ( Ui)y ,  over  Ahex  would  cancel  out  and  be  null.  The 
same  feature  has  been  used  in  the  development  of  the  original  second-order  CESE  method.  Therefore,  for 
the  fourth-order  CESE  method,  the  space-time  flux  passing  through  the  top  surface  of  the  CE,  i.e.,  Eqs.  (20) 
and  (21),  is  an  implicit  function  of  (ui) 3x  and  its  second-  and  third-order  derivatives.  However,  as  will  be 
shown  in  the  following,  these  unknowns  do  not  need  to  be  solved  simultaneously.  Instead,  we  will  solve  for 
the  second-  and  third-order  derivatives  of  iq  first.  Then,  Eq.  (21)  becomes  an  explicit  function  of  (iq) 3x 
only. 

To  proceed,  we  integrate  three  additional  model  equations  shown  in  Eq.  (13)  to  calculate  (ui)xx,  ( Ui)xy , 
and  (ui)yy  and  their  first-order  derivatives,  i.e.,  the  third-order  derivatives  of  ut ,  at  point  3X.  To  this  end, 
we  simply  apply  the  original  second-order  CESE  method  three  times  for  space-time  flux  conservation  of 
{ui)xx,{ui)Xy,  and  (itj) yy.  Although  flux  conservation  of  these  three  variables  are  not  prescribed  by  the 
original  model  equations,  the  invoked  assumption  of  smoothness  and  differentiable  of  the  primary  unknowns 
Ui  up  to  the  third-order  would  ensure  the  validity  of  these  three  additional  conservation  equations.  The 
procedure  of  solving  the  second-  and  third-order  derivatives  of  Ui  is  identical  to  that  of  the  original,  second- 
order  CESE  method  for  solving  m  and  its  first-order  derivatives.  For  conciseness,  no  discussion  about  the 
integration  procedure  is  repeated  here. 

With  the  known  second-  and  third-order  derivatives  of  (ui)3x,  Eq.  (21)  becomes  an  explicit  function  of 
the  conserved  variable  of  ( Ui)3x  at  the  solution  point, 

(ffi) 3X  =  ~7  [(-Fi)side  T  (Fi) bot  T  (-Fi)top,2nd,3rd]  (22) 

-^hex 

where  (Fi)top,2nd,3rd  denotes  the  part  of  the  space-time  flux  in  Eq.  (21),  which  can  be  calculated  by  using 
the  known  second-  and  third-order  derivatives  of  rq. 

The  integrals  in  Eqs.  (15),  (18),  and  (20)  can  be  explicitly  calculated  in  several  different  manners. 
In  this  work,  we  calculate  the  integrals  by  using  a  symbolic  mathematic  package.  As  such  all  unknowns 
in  the  fourth-order  CESE  method  are  calculated  except  the  first-order  derivative  of  Uj.  The  algorithm 
of  calculating  the  first-order  derivatives  is  identical  to  that  of  the  second-order  CESE  method.  Essentially, 
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based  on  the  known  even-order  derivatives  (one  order  lower  than  the  desired  odd-order  derivatives),  a  central- 
differencing  procedure  is  applied  to  the  even-order  derivatives  to  calculate  the  odd-order  derivatives.  The 
same  algorithm  is  used  to  calculate  third-order  derivatives  based  on  the  known  second-order  derivatives  in 
the  above  fourth-order  CESE  method.  For  completeness,  the  central  difference  method  is  illustrated  in  the 
following  subsection. 


II.  D. 4-  Odd- Order  Derivatives 

The  odd-order  derivatives  of  the  conserved  variables  are  calculated  by  using  a  central-differencing  method, 
which  is  identical  to  that  used  in  the  original,  second-order  CESE  method.  (1  The  mesh  stencil  employed  for 
calculating  the  odd  derivatives  is  shown  in  Fig.  5,  in  which  symbol  squares  denote  vertices  of  triangles,  crosses 
denote  solution  points,  and  plus  signs  the  shifted  locations  of  the  solutions  points  for  the  CFL  insensitive 
schemes.  The  spatial  location  of  shifted  solution  points  (the  plus  sign)  are  calculated  by 

(xj+,yj+)  =  [r(a;3x  -  Xjx )  +  xjx ,  r(y3x  -yjx)  +  yjx],  (23) 


where  j  =  1,  5,  7  indicate  the  three  neighbors  of  the  CE  centered  at  point  3X  where  the  solution  at  a  new 
time  step  is  sought,  r  is  a  function  of  the  local  CFL  number,  (Xjx,yjx),  with  j  =  1,5,7,  are  the  original 
solution  points  of  the  three  neighbors,  and  (xj+ ,  yj+),  with  j  =  1,5,7,  are  the  shifted  solution  points  of  the 
three  neighbors.  The  solutions  at  the  previous  time  step  n  —  1/2  to  be  used  for  the  time  marching  calculation 
are  stored  at  ( Xjx  ,  y,j x ),  with  j  =  1,  5,  7.  Thus,  we  first  take  a  Taylor  expansion  from  (x? x  ,  y^x  )  at  the  time 
level  n  —  1/2  to  the  shifted  locations  (Xj+ ,  y]+ )  at  the  new  time  level  n: 


A  A— a A—a—b 


[(«i) 


xTyJ  \jJ, 


EE 

a— 0  6=0 


dBi 


E'  \dxI+adyJ+bdtk+c 


alblcl 


(xJ+  -  xjx)a  (yj+  -  y0x)°  (  .  (24) 


Aided  by  Eq.  (24)  solutions  of  unknowns  at  the  three  shifted  solution  points  in  the  new  time  level  are 
calculated.  The  four  points  1+,  5+,  7+,  and  3X  form  three  triangles.  Consider  triangle  A3xl+5+.  The 
spatial  derivatives  of  the  unknowns  in  this  triangle  can  be  readily  calculated: 

xi+  ~  x3*  Vi+  ~  2/3x 
x5+  ~  X3x  2/5+  —  2/3x 


(ui)  a 

(u%)x 


[K) 


xIyJ\ i+ 


-  [K)2 


x'y-'  J  3X 


(25) 


Similar  calculation  as  that  in  Eqs.  (25)  and  (24)  can  be  performed  to  find  the  odd-order  derivatives  [(ui)xi+iyj , 
(ui)xiyj+ 1]  in  triangles  A3x7+5+  and  A3xl+7+.  For  shock  capturing,  these  three  sets  of  odd-order  deriva¬ 
tives  are  weighted  by  using  the  following  formula: 


( d‘i)xI+1yJ 
(^+)  x1  yt7+1 


=  w 
=  w 


{[{Ui)xI+1yJ 

r. 

[(ui)Xi+1yJ]^  ^ 

[(' Ui)xi+iyj]^  ,aj 

(26) 

r, 

(27) 

where  the  superscript  (r)  with  r  =  1,2,3  denotes  the  three  triangles  A3xl+5+,  A3x7+5+a,  and  A3xl+7+, 
and  the  weighting  function  W  is  defined  as 


W  = 


(6*26*3)a/z  +  (6?,6i)acj  +  ( 9i02)ail) 

(pxo2)a  +  {e2e3)a  +  (030i)a 


where  8  is  the  magnitude  of  the  odd-order  derivative  vector: 


9r 


1  (0 


J  +{[(«i)a 


(28) 


(29) 


Equations  (24-??)  are  used  to  calculate  the  odd-order  derivatives  of  the  conserved  variables  to  the  new 
time  step.  Shown  in  the  above  equations,  the  calculation  of  the  odd-order  derivatives  depend  on  the  even- 
order  derivatives  of  the  next  lower  order.  To  recapitulate,  the  algorithm  of  the  fourth-order  CESE  method 
can  be  organized  in  the  following  steps:  (i)  Apply  the  original  second-order  CESE  method  three  times  to 
solve  for  the  additional  equations  in  terms  of  the  second-order  derivatives  of  Ui  as  shown  in  Eq.  (13)  for 
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the  solution  of  the  second-order  derivatives  based  on  the  space-time  flux  conservation,  (ii)  Calculate  all 
the  third-order  derivatives  by  central  differencing  the  solutions  of  the  second-order  derivatives,  (iii)  Apply 
the  newly  developed  fourth-oder  CESE  method  to  calculate  the  zero-order  derivatives,  or  the  conserved 
variables  themselves,  (iv)  Finally,  we  calculate  the  first-order  derivatives  by  central  differencing  the  conserved 
variables. 

The  above  recursive  algorithm  can  be  extended  to  sixth-  or  even  higher-order  methods.  The  mesh 
stencils  employed  for  all  high-order  methods  is  identical  to  that  used  for  the  original  CESE  method.  For 
the  sixth-order  CESE  method,  we  first  repeatedly  apply  the  second-order  CESE  method  to  solve  for  all  the 
fourth-order  derivatives  of  14.  Then,  the  fourth-orders  derivatives  are  central  differenced  to  obtain  the  fifth- 
order  derivatives.  Next,  we  would  apply  the  fourth-order  method  three  times  to  solve  for  all  the  second-order 
derivatives,  followed  by  central  difference  procedure  for  the  third-order  derivatives.  Finally,  one  can  apply 
the  six-order  CESE  method  to  calculate  the  conserved  variables,  followed  by  central  difference  procedure  for 
the  second-order  derivatives. 


III.  Results  and  Discussion 


In  this  section,  we  report  the  numerical  results  obtained  by  applying  the  newly  developed,  two-dimensional, 
fourth-order  CESE  method.  Three  sets  of  model  equations  are  considered,  including  (i)  a  scalar  convection 
equation,  (ii)  the  linearized  Euler  equations  for  acoustics,  and  (iii)  the  Euler  equations.  Numerical  results 
of  four  cases  are  reported:  one  case  each  for  solving  the  scalar  equation  and  the  acoustic  equations,  and  two 
cases  for  solving  the  Euler  equations.  In  all  four  cases,  the  third-order  Taylor  series  is  used  to  discretize  the 
conserved  unknowns  and  the  fluxes  inside  each  SE.  Thus,  the  numerical  scheme  employed  is  fourth-order  in 
accuracy.  To  determine  the  convergence  rates,  we  measure  the  L2  norm  of  errors.  For  a  conserved  variable 
Ui,  the  nth  norm  Ln  with  n  =  1  or  2,  is  defined  as 


L 


n  — 


/5Z  [(U*)l>num 


(wj)j',ana]  Ay, 


where  the  subscript  num  denotes  the  numerical  result,  ana  the  analytical  solution,  and  Aj  is  the  area  of  cell 
j.  The  order  of  accuracy  of  the  numerical  results  is  assessed  by  the  convergence  test.  We  repeat  the  same 
simulation  with  several  sets  of  meshes  with  decreasing  cell  sizes.  The  convergence  tests  are  performed  for 
the  first  three  cases,  which  involve  linear  solutions.  In  the  fourth  case,  the  solution  involves  shock  waves.  In 
the  vicinity  of  shocks,  a  reweighing  function  is  activated  and  artificial  damping  is  added  to  solutions  so  that 
oscillations  of  solution  in  the  near  shock  region  would  be  suppressed.  Thus  the  order  of  accuracy  would  be 
reduced  to  be  first-order  around  the  shock  wave. 

Each  convergence  test  is  done  by  using  two  types  of  meshes:  (i)  uniform  meshes  that  is  created  by  splitting 
a  quad  mesh  into  4  triangles,  (ii)  unstructured  meshes  composed  of  random  triangles.  We  used  Cubit3,  a 
DoE  package  developed  by  using  the  advancing-front  algorithm,  to  generate  all  meshes.  For  the  uniform 
meshes  the  option  of  QTri  mesh,  provided  by  Cubic  was  used.  For  the  unstructured  meshes  the  option  of 
TriAdvance  was  used.  In  using  the  TriAdvance  algorithm,  we  first  employed  the  Delaunay  algorithm  to 
generate  the  preliminary  meshes.  We  then  use  the  QTri  method  to  refine  the  meshes.  The  remainder  of  this 
section  reports  the  numerical  results  of  the  four  cases. 


III. A.  The  Scalar  Advection  Equation 

The  equation  of  concern  is 


^  +  ax-+a,--°. 

A  square  domain  —  tt  <  x  <  n,  —7 r  <  y  <  it  is  considered.  Periodic  boundary  conditions  are  imposed  to 
all  four  sides  of  the  computational  domain.  The  initial  condition  is  u{x,  y,0)  =  sin(axx  +  avy).  Due  to 

the  periodic  condition,  the  analytical  solution  ua na  =  sin(axx  +  avy  +  att ),  with  at  =  +  ay-  In  the 

present  calculation,  we  let  ax  =  ay  =  1.  In  calculation,  we  let  the  simulation  run  until  t  =  15,  and  the  CFL 

ahttp:  / /cubit. sandia.gov/index.html 


9  of  20 


American  Institute  of  Aeronautics  and  Astronautics 


number  «  1.  Table  2  shows  the  calculated  convergence  rate  by  using  the  newly  developed  fourth-order  CESE 
method.  Table  2(a)  shows  the  results  of  using  uniform  meshes  and  the  convergence  up  to  fifth  order  is  found. 
Table  2(b)  shows  results  of  using  unstructured  triangular  meshes  and  the  convergence  rate  is  fourth-order. 
When  calculating  the  convergence  rates  for  the  uniform  meshes,  the  characteristic  length  employed  is  the 
square  root  of  the  area  of  each  triangular  cell.  For  unstructured  meshes,  it  is  the  square  root  of  the  area  of 
the  largest  triangular  cell  in  the  mesh. 

Table  2:  Convergence  rates  of  the  fourth-order  CESE  method  for  solving  the  scalar  advection  equation. 

(a)  Structured  grid 


number  of  cells 

h 

L2  error 

O  (L2) 

400 

9.870E-02 

2.158E-01 

- 

1600 

2.467E-02 

6.726E-03 

5.00 

3600 

1.097E-02 

8.762E-04 

5.03 

6400 

6.169E-03 

2.065E-04 

5.02 

10000 

3.948E-03 

6.735E-05 

5.02 

(b)  Unstructured  grid 

number  of  cells 

h 

L2  error 

O  (L2) 

232 

1.690E-01 

2.929E-01 

- 

926 

4.255E-02 

9.727E-03 

4.94 

2068 

1.909E-02 

1.315E-03 

4.99 

3694 

1.068E-02 

3.131E-04 

4.94 

5774 

6.835E-03 

1.022E-04 

5.01 

III.B.  The  Linear  Acoustic  Equations 


The  acoustic  equations  are  the  linearized  Euler  equations  for  compressible  flows: 


dp 

dt 


Po 


dvx 

dx 

dvx 

dt 


dvy_ 

dt 


dvx\ 
dx  J 


=  0 


aod£_  =  Q 
po  dx 

Po  dy 


where  p,  vx,  vy,  a o,  and  po  are  respectively  the  density,  x  and  y  component  of  the  velocity,  the  speed  of 
sound  of  the  free  stream,  and  the  density  of  the  free  stream.  A  square  domain  is  considered:  —2tt<x<2tt 
and  —2ir  <  y  <  2n.  We  then  impose  the  periodic  boundary  conditions  to  the  four  sides  of  the  computational 
domain.  When  t  =  0,  the  initial  conditions  are 


p(x,  y,  0)  =  po  +  p'  sin(axx  +  avy) 

vx(x ,  y,  0)  =  -Oq  —  —  sin(axx  +  avy) 
po  at 

vy(x,y,  0)  =  -al  —  —  sin(axx  +  avy), 

Po  at 

where  p'  is  a  prescribed  amplitude  of  the  perturbation.  As  such,  the  analytical  solution  is 


Pana  =  Po  +  p'  sin(axx  +  dyy  +  att) 

(w)ana  =  -a%  —  —  sin(axx  +  dyy  -I-  att) 
Po  at 

(^y)ana  Uq  —  sin(uxx  T  ayy  att), 
Po  at 
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where  at  =  —ao^Ja%.  +  a y.  In  this  test  case,  we  let  ax  =  ay  =  1,  p'  =  0.01,  ao  =  Po  =  1-0.  The  time-marching 
calculation  is  run  until  t  =  18.  The  CFL  number  is  about  0.7  for  the  case  using  the  unstructured  mesh,  and 
about  0.75  for  the  case  using  the  uniform  mesh. 

Shown  in  Fig.  6,  the  convergence  rate  is  approximately  4.3  for  the  unstructured-mesh  case,  and  4.6 
for  the  uniform-mesh  case.  The  convergence  rate  for  the  non-uniform  case  includes  the  four  most  courses 
simulations.  In  the  figure,  h  denote  the  length  scale  of  a  typical  mesh  cell,  i.e. ,  h  ss  Ax  ~  Ay.  For  similar 

h,  the  result  obtained  by  using  the  unstructured  mesh  tends  to  be  more  accurate  than  that  by  using  the 
structured  mesh.  This  is  true  up  until  an  h  value  of  0.1.  At  this  point  the  L2  norm  for  the  unstructured 
mesh  begins  to  increase.  The  sudden  decrease  in  the  convergence  rate  was  further  investigated  by  the  Euler 
equation. 

III.C.  The  Euler  Equations  for  Linear  Waves 

This  section  reports  numerical  results  of  the  two-dimensional  Euler  equations  obtained  by  using  the  fourth- 
order  CESE  method.  Two  cases  are  included.  The  first  one  is  related  to  linear  waves,  which  are  used  to 
determine  the  convergence  rate  of  the  solver.  The  second  test  shows  the  shock  capturing  capabilities  of  the 
new  solver.  For  the  linear  wave  case,  the  following  initial  conditions  are  used: 

p{x,  y,  0)  =  1.0  +  0.002  sin(axx  +  avy  +  aR) 
vx(x,y,  0)  =  0.5  vy(x,y,0)  =  0.5  p(x,  y,  0)  =  24/3/y 

where  ax  =  ay  =  n  and  at  =  —  {ax Vx  +  ayVy).  The  computational  domain  is  2  <  x  <  2,  —2  <  y  <  2.  The 
calculation  is  done  for  time  period  of  0  <  t  <  4.  The  total  time  period  of  4  corresponds  to  2  periods  of  wave. 
The  periodic  boundary  condition  is  imposed  to  the  four  sides  of  the  square  domain. 

For  the  convergence  test,  the  numerical  solution  of  the  fourth-order  CESE  method  was  compared  with 
that  of  the  second-order  CESE  method.  Shown  in  Fig.  7,  the  second-order  simulation  achieves  a  convergence 
rate  of  3.1,  while  the  fourth-order  simulations  achieves  a  convergence  rate  of  4.1.  When  the  Li  norm  is 
used,  the  convergence  rates  decrease  by  one  for  both  cases.  When  the  Loo  norm  is  used,  the  convergence 
rate  increases  by  one.  For  all  convergence  rates,  the  characteristic  size  h  is  equal  to  the  square  root  of  the 
averaged  areas  of  all  cells.  In  all  cases,  the  CFL  number  varies  between  0.5  and  0.8. 

III.D.  The  Euler  Equations  for  an  Expanding  Shock 

In  this  case,  the  new  fourth-order  CESE  method  has  been  used  to  simulate  a  compressible  flow  with  shock 
wave.  The  computational  domain  is  a  square.  The  length  of  each  side  is  2.  Non-reflecting  boundary  condition 
is  imposed  to  all  four  sides  of  the  computational  domain.  The  initial  conditions  is  a  two-dimensional  set  up 
of  a  shock  tube.  Inside  the  radius  r  <  Rq,  pressure  and  density  are  higher  than  that  outside  of  the  radius, 

i. e.,  r  >  Rq. 


_  I  1.0, 1.0  if  r  <  R0 
P,P  ~  jo. 125, 0.1  if  r  >  Rq’ 

where  r  =  y/(x  —  Xq)2  +  (y  —  yo)2,  Ro  =  0.5,  and  (xo,  yo)  is  the  center  of  the  high-pressure  domain.  At  the 
beginning  of  the  simulation,  the  diaphragm  separating  the  high  pressure  region  and  the  low  pressure  region 
is  lifted  and  a  shock  wave  expands  outwards  and  a  expansion  wave  converges  towards  the  center.  This  case 
was  simulated  by  using  both  the  second-  and  fourth-order  CESE  solvers.  When  running  the  fourth-order 
simulations  it  was  found  that  the  current  suite  of  limiters  commonly  used  in  the  second-order  CESE  method 
were  not  potent  enough  to  maintain  stability. 

A  new  limiter  was  developed  that  insures  that  the  Taylor  series  convergences.  This  is  accomplished  by 
looking  at  the  influence  that  the  next  higher  derivatives  has  on  the  summation  of  the  lower  derivatives  in 
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the  Taylor  series. 


/  Qui+J+ i  \ 
\dxI+1dyJ  ) 


hi+j 

i\j\ 


and 


/  QuI+J+1  \ 
\dxIdyJ+1 ) 


h<P 


hi+j 

i\j\ 


where  h  is  the  square  root  of  the  cell  area,  /?  is  an  adjustable  parameter  on  the  order  of  1.  If  this  check  does 
not  pass,  the  higher-order  derivative  is  set  to  be  null.  This  check  is  done  twice.  First  after  the  second  and 
third  derivatives  are  calculated,  then  again  after  the  zero  and  first  derivatives  are  calculated.  It  should  be 
noted  that  this  test  is  done  in  addition  to  the  weighted  odd  derivative  calculation  expressed  in  Eq.  28. 

This  case  of  shock  wave  was  tested  by  using  three  different  meshes  with  5770,  23000  and  91754  cells.  The 
results  can  be  seen  in  Fig.  10.  The  results  from  the  fourth-order  CESE  scheme  were  compared  to  the 
second-order  CESE  scheme.  The  plots  shown  in  Fig.  10  are  the  values  plotted  along  a  line  going  from  the 
lower  left  to  the  upper  right  corners.  In  all  cases  the  fourth-order  solution  has  more  dissipation  near  the 
middle  of  the  domain.  For  the  5770  cell  case  the  second-order  version  was  not  able  to  produce  the  correct 
shock  and  contact  location,  while  the  fourth-order  version  was  able  to.  For  the  other  two  cases  the  results 
are  very  similar  with,  the  fourth-order  version  has  slightly  sharper  contacts  but  more  dissipation  towards 
the  center  of  the  domain.  This  dissipation  decreases  as  the  resolution  increases.  In  all  cases  the  value  of  /? 
was  1  and  a  was  2.  The  average  CFL  number  were  0.84,  0.88  and  0.93  for  the  5770,  23000  and  91754  cases 
respectively. 


IV.  Concluding  Remarks 

This  paper  reports  a  new  fourth-order  CESE  method  for  solving  the  two-dimensional  linear  and  nonlinear 
hyperbolic  partial  differential  equations  by  using  unstructured  meshes.  The  method  is  an  extension  of 
Chang’s  high-order  CESE  methods  for  one-dimensional  equations.  The  method  here  focuses  on  the  use 
of  unstructured  triangular  meshes  for  solving  a  set  of  coupled  equations.  The  new  fourth-order  CESE 
method  retains  all  favorable  features  of  the  original  second-order  CESE  method,  including  (i)  the  use  of 
the  most  compact  mesh  stencil  involving  only  the  immediate  neighboring  nodes  of  the  central  node  where 
the  unknowns  are  sought,  (ii)  the  stability  constraint  of  the  four-order  CESE  method  remains  to  be  CFL 
<  1,  and  (iii)  completely  explicit  operation  in  the  time  marching  calculation.  Moreover,  the  presented 
formulation  is  general  and  recursive;  it  can  be  easily  extended  to  the  sixth-  or  a  even  higher-oder  CESE 
method.  To  demonstrate  the  capabilities  of  the  new  method,  we  have  reported  numerical  solutions  of  three 
different  sets  of  model  equations,  including  the  scalar  advection  equation,  the  acoustic  equations,  and  the 
Euler  equations.  For  solutions  of  linear  problems,  we  demonstrated  fourth-order  convergence  rates  by  using 
completely  unstructured  meshes.  When  using  a  uniform  mesh,  numerical  results  show  super  convergence. 
We  have  also  solved  the  Euler  equations  for  modeling  supersonic  flows  with  shocks.  Numerical  results  of 
an  expanding  shock  demonstrated  that,  by  using  the  original  reweighing  function  as  that  in  the  original 
second-order  CESE  method  to  the  first-  and  third-order  derivatives,  the  new  fourth-order  CESE  method  can 
accurately  capture  moving  shock  waves. 
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8  2  9 


10 


Fig.  1:  The  CESE  domain  using  a  triangular  mesh  with  points  of  interest  marked  by  squares(cell  vertices), 
circles(cell  centroid)  and  crosses(CE  centroid/solution  point). 


(a)  Conservation  Element 


(b)  Solution  Element 


Fig.  2:  The  CE  and  SE  centered  on  the  centroid  mesh  point  (j,n). 
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V.  Figures 


(X2,y2,tn)r 


>2,2/2,*"  1/2)'r 


Fig.  3:  The  side  faces  associated  with  a  solution  point 
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o.i 


non-uniform  +  uniform 

Fig.  6:  Convergence  test  for  the  linear-acoustic  equations. 
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h 


2nd  4  4th  x 

Fig.  7:  Convergence  test  for  the  Euler  Equations 
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Fig.  8:  Fourth  order  convergence  test  for  the  Euler  Equations  on  an  unstructured  mesh. 
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Fig.  9:  An  examle  domain  used  in  the  shock  case 
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Fig.  10:  Comparison  of  the  fourth-  and  second-order  CESE  schemes  for  capturing  shocks  at  three  different 
resolutions. 
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